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Abstract — The  emerging  memristor  devices  have  recently  re¬ 
ceived  increased  attention  since  HP  Lab  reported  the  first  TiCL- 
based  memristive  structure.  As  it  is  at  nano-scale  geometry  size, 
the  uniformity  of  memristor  device  is  difficult  to  control  due 
to  the  process  variations  in  the  fabrication  process.  The  in¬ 
curred  design  concerns  in  a  memristor-based  computing  system, 
e.g,  neuromorphic  computing,  can  be  very  severe  because  the 
analog  states  of  memristors  are  heavily  utilized.  Therefore,  the 
understanding  and  quantitative  characterization  of  the  impact 
of  process  variations  on  the  electrical  properties  of  memristors 
become  crucial  for  the  corresponding  VLSI  designs.  In  this  work, 
we  examined  the  theoretical  model  of  TiCL  thin-film  memristors 
and  studied  the  relationships  between  the  electrical  parameters 
and  the  process  variations  of  the  devices.  A  statistical  model 
based  on  a  process-variation  aware  memristor  device  structure 
is  extracted  accordingly.  Simulations  show  that  our  proposed 
model  is  3  ~  4  magnitude  faster  than  the  existing  Monte- 
Carlo  simulation  method,  with  only  ~  2%  accuracy  degradation. 
A  variable  gain  amplifier  (VGA)  is  used  as  the  case  study  to 
demonstrate  the  applications  of  our  model  in  memristor-based 
circuit  designs. 


I.  Introduction 

Based  on  the  completeness  of  circuit  theory,  Prof.  Chua 
predicted  the  existence  of  the  fourth  fundamental  passive 
circuit  element  -  memristor  in  1971  [1].  Till  2008,  a  real 
memristor  device  was  implemented  through  a  Ti02  thin-film 
structure  [2].  Since  then,  many  different  types  of  memristor 
materials  and  devices  have  been  reported  or  rediscovered. 

Memristors  show  many  promising  characteristics  as  the 
next-generation  data  storage  devices,  such  as  non- volatility, 
low-power  consumption,  high  integration  density  and  excellent 
scalability  [3]  [4].  Also,  the  special  property  that  remembers 
historical  profile  of  the  electrical  excitations  on  the  device  [5] 
makes  memristor  an  ideal  candidate  to  realize  the  synapse 
behavior  in  electronic  neural  networks  [6]  [7]. 

As  process  technology  scales  down  to  45nm  and  below,  pro¬ 
cess  variations,  e.g.,  line-edge  roughness  (LER),  oxide  thick¬ 
ness  fluctuation  (OTF),  and  random  discrete  doping  (RDD)  [8], 
lead  to  significant  device  parameter  fluctuations.  The  impact  of 
process  variations  on  a  memristive  system  may  be  more  severe 
than  a  conventional  digital  design  because  of  the  utilization 
of  the  analog  states  of  memristors.  Many  models  have  been 
carried  out  for  various  memristor  devices,  including  Ti02  thin- 
film  memristors  [9],  spintronic  memristors  [10],  and  ion  con¬ 


ductor  chalcogenide-based  memristors  [11].  Very  recently,  Hu 
et  al.  presented  a  geometry-variation  aware  memristor  device 
structure  model  [12].  The  associated  Monte-Carlo  simulation 
flow,  however,  is  very  time-consuming. 

In  this  work,  we  developed  a  fast  statistical  model  to 
simulate  the  electrical  properties  of  Ti02  thin-film  memris¬ 
tors.  Starting  with  the  theoretical  model  of  a  Ti02  thin-film 
memristor,  we  explored  the  influences  of  geometry  variations 
on  the  electrical  parameters  of  the  device.  On  top  of  that,  a 
statistical  model  with  the  merits  of  both  high  accuracy  and 
low  runtime  cost  is  proposed.  Compared  to  the  existing  work 
in  [12],  our  model  significantly  improves  the  runtime  cost  by 
3^4  orders  of  magnitude;  and  reduces  the  input  data  set 
down  to  a  few  variables.  The  simulation  accuracy  maintains 
within  ~  2%  discrepancy  in  the  whole  working  region  of  the 
memristor  device,  compared  to  the  results  using  the  model 
in  [12].  Furthermore,  we  demonstrated  a  design  example  -  an 
op-amp  based  variable  gain  amplifier  (VGA),  and  analyzed  the 
impact  of  process  variations  to  the  design  with  the  extracted 
statistical  model. 

The  rest  of  paper  is  organized  as  follows:  Section  II  in¬ 
troduces  the  memristor  basics,  including  the  theory  and  the 
physical  mechanisms  of  Ti02  thin-film  memristors,  and  ana¬ 
lyzes  the  impact  of  process  variations  on  the  device  electrical 
properties;  Section  III  illustrates  the  theoretical  fundamentals 
of  our  proposed  model;  Section  IV  explains  the  methodology 
to  generate  our  proposed  statistical  model;  Section  V  demon¬ 
strates  the  effectiveness  and  accuracy  of  our  proposed  model 
by  some  examples;  and  Section  VI  concludes  our  work. 


II.  Background 


A.  Memristor  Theory 


The  original  definition  of  the  memristor  is  derived  from  the 
completeness  of  circuit  theory:  besides  resistor,  capacitor  and 
inductor,  there  exists  the  fourth  basic  two-terminal  element  that 
uniquely  defines  the  relationship  between  the  magnetic  flux 
(</?)  and  the  electric  charge  ( q )  passing  through  the  device, 
or  df  =  M  •  dq  [1].  Since  f  =  f  Vdt  and  q  =  f  Idt ,  the 
definition  of  the  memristor  can  be  generalized  as: 


V  =  M(u,I)-I 

f  =  /(“>,/) 


(i) 
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C.  Process -Variation  Aware  Monte-Carlo  Simulations 


where  uj  is  a  state  variable;  M(cc,  I)  denotes  the  instantaneous 
memristance,  which  may  vary  over  time.  Under  certain  condi¬ 
tions,  M(cu,  I)  and  /(cc,  I)  may  be  simplified  as  the  explicit 
functions  of  I. 


B.  Ti02  Memristor 

In  2008,  HP  Lab  demonstrated  the  memristive  effect  by 
moving  the  doping  front  along  a  Ti02  thin-film  [2],  Fig. 
1(a)  illustrates  the  Pt/TiC^/Pt  structure.  The  stoichiometric 
TiCU  with  exactly  2:1  ratio  of  oxygen  to  titanium  has  a  low 
conductivity  and  can  be  considered  as  an  insulator.  When 
the  Ti02  loses  a  certain  amount  of  oxygen  (called  oxygen- 
deficient  titanium  dioxide),  its  conductivity  becomes  relatively 
high  as  a  semiconductor.  A  positive  bias  causes  the  oxygen 
vacancies  to  drift  to  the  pure  Ti02  region  and  lowers  the 
overall  resistance  continuously.  A  negative  bias,  however, 
reverses  the  above  process  and  raises  the  overall  resistance 
[3]. 

Fig.  1(b)  illustrates  the  equivalent  circuit  of  such  a  TiC>2 
thin-film  memristor  as  two  resistors  connected  in  series.  Rl 
( Rh )  denotes  the  lowest  (highest)  resistance  of  the  device 
when  the  memristor  is  fully  doped  or  undoped,  respectively. 
The  total  memristance  of  a  TiCU  memristor  can  be  expressed 
as 

M(a)  =  Rl  •  cx  +  Rh  •  (1  -  cy).  (2) 

Here  a  (0  <  a  <  1)  is  the  relative  doping  front  position,  which 
is  the  ratio  of  doping  front  position  over  the  total  thickness  of 
the  Ti02  thin-film. 

The  velocity  of  doping  front  movement  v{t ),  which  is  driven 
by  the  voltage  applied  across  the  memristor  V{t)  can  be 
expressed  as 

V(t)  =  ^  =  P'v  •  •  M(o!)  ’  (3) 

where  pv  is  the  equivalent  mobility  of  dopants;  h  is  the  thick¬ 
ness  of  the  TiC>2  thin  film;  and  M(a)  is  the  total  memristance 
when  the  relative  doping  front  position  is  a.  The  change  of  a 
over  time  t  can  be  obtained  by  solving  the  differential  equation 
Eq.  (3)  as 

a(t)  =  (4) 

Here  A  =  pv  •  ^  •  f^V  ( t)dt  exhibits  the  memristor’ s  historic 
behavior;  and  B  =  RH  •  a0  +  §  •  (Rl  —  Rh)  •  *s  determined 
by  the  initial  condition  cxq  =  a(to). 


The  process  variations  significantly  influence  the  electrical 
properties  of  nano-devices.  For  a  TiC>2  thin-film  memristor, 
there  are  three  major  sources  of  process  variations:  (1)  Line 
edge  roughness  (LER),  which  is  generated  in  the  lithography 
and  etching  steps  and  causes  the  deviation  of  the  printed 
image’s  edge  from  its  designed  pattern;  (2)  Oxide  thickness 
fluctuation  (OTF),  which  is  generated  in  the  deposition  pro¬ 
cess;  and  (3)  Random  discrete  doping  (RDD),  which  leads 
to  the  material  non-uniformity  within  the  device.  In  the  rest 
of  paper,  without  loss  of  generality,  we  use  c o  and  uJ  to 
respectively  represent  the  design  value  and  the  actual  value 
under  process  variations  for  any  given  variable  c o. 

Very  recently,  a  Monte-Carlo  simulation  method  was  pro¬ 
posed  to  analyze  the  impact  of  geometry  variations  on  the 
electrical  properties  of  TiC>2  thin-film  memristors  [12].  A  3D 
device  structure  including  the  geometry  variation  information 
is  generated  by  performing  a  sanity  check  of  the  device 
characterization  parameters.  In  the  Monte-Carlo  simulations, 
the  memristor  device  is  divided  into  many  small  filaments 
sandwiched  between  two  electrodes.  Within  a  filament  i,  the 
cross-section  area  s'  =  I'j  •  z[  and  thickness  h[,  can  be 
considered  as  constants,  whose  value  can  be  modeled  by 
taking  into  account  the  effects  of  LER  or  OTF,  respectively. 
The  corresponding  memristance  and  doping  front  position  of 
filament  i  can  be  represented  as 


■)  =  R'L,i  ■  ai  +  RH,i  ■  (!  -  «i). 


(5) 


and 


a'M 


R'h  i—R'r 


(6) 


r'  t 

where  the  coefficients  A[  =  /jlv  •  •  fto  V  ( t )  •  dt  and  B[  = 


R'ha  *  ao  +  2  (R'la  R'h  a) 

'  aQ. 

The  two  boundary  resistances  of  filament  i 

•  r'l,i  anci  R'n.p 

which  can  be  calculated  by 

RL,i  =  fo 

ki  Pon  . 

1  s'(a')  aOLV 

(7) 

and 

R'h,i  ~  fc 

hi  Poff  .  Anl 
)  sj(a')  aaR 

(8) 

vary  from  filament  to  filament  under  the  process  variations. 
Here  pon  and  pQff  are  the  electrical  resistivity  of  the  doped 
and  updoped  Ti02  thin-film,  respectively,  which  suffer  from 
the  RDD  issue. 

The  overall  instantaneous  memristance  of  the  Ti02  memris¬ 
tor  can  be  calculated  as  the  total  resistance  of  all  n  filaments 
connected  in  parallel  as 


M'(t)  —  Sn=i\ /M/-  (9) 

Due  to  the  thin  film  thickness  fluctuation  (and/or  the  rough¬ 
ness  of  the  electrode  contact),  the  doping  front  movement 
inside  every  filament  varies.  Therefore,  a  concept  named 
“average  doping  front  position”  is  introduced  as 


a'{t) 


(10) 


2 

346 


- M-C  simulation - Norm-fit 


In  the  Monte-Carlo  simulation  flow  proposed  in  [12],  a  fil¬ 
ament  based  process  variation  model  is  carried  out.  However, 
due  to  its  large  runtime  cost,  this  technique  is  infeasible  in 
large  scale  circuit  and  system  designs.  In  this  work,  we  will 
proposed  a  more  efficient  process-variation  aware  memristor 
model  with  a  very  slight  accuracy  relaxation. 

III.  An  Analysis  Of  Memristor  Model 

Compared  to  other  fundamental  passive  devices,  the  model¬ 
ing  of  memristors  has  some  unique  challenges:  First,  as  a  time- 
varying  device  whose  properties  are  determined  by  the  historic 
profile  of  electrical  excitations,  not  only  the  variations  of  the 
start  and  the  end  states,  but  also  that  of  the  intermediate  time- 
varying  state  must  be  modeled;  Second,  as  an  analog  device 
that  can  be  used  in  diverse  circuit  styles  and  environments, 
various  device  properties,  e.g.,  the  total  memristance  and  its 
changing  rate,  the  timing  and  frequency  responses,  and  I-V 
characteristic,  must  be  modeled;  finally,  the  calibration  and 
validation  of  the  proposed  models  can  be  challenging  due  to 
the  lack  of  published  measurement  data. 

When  the  device  variations  are  within  the  reasonable  range, 
we  can  assume  the  ratio  between  the  actual  device  parameter 
and  its  designed  value  as  a  polynomial  expression,  or  x'  = 
r]-x.  Here  77  is  a  coefficient  representing  the  effects  of  process 
variations.  Our  goal  is  to  find  an  efficient  methodology  to 
compute  the  variation-aware  coefficient  77. 

The  total  memristance  M',  which  is  a  time- varying  pa¬ 
rameter,  can  be  uniquely  defined  by  R'H ,  R'L ,  and  In 

this  section,  we  will  examine  the  variations  of  R'H,  R'L ,  and 
a'  (t)  first,  and  then  derive  the  corresponding  process-variation 
aware  memristance  model.  Table  I  summarizes  the  designed 
values  of  the  Ti02  memristor  geometry  dimensions  and  its 
electrical  parameters  adopted  in  our  work  [2]. 

A.  Distribution  of  R'H  and  R'L 

For  a  given  Ti02  memristor,  we  use  R'H  and  R'L  to  denote 
its  actual  highest  and  lowest  total  memristances,  respectively. 
As  Eq.  (7)  and  (8),  the  geometry  variations  {h[  and  s')  influ¬ 
ence  both  R'h  ■  and  R'L  ■  simultaneously  within  the  filament 
i.  It  indicates  R'H  ■  and  R'L  ■  are  correlated.  However,  they  are 
not  fully  correlated  because  of  the  randomness  in  pQff  and 
pon,  which  is  incurred  by  RDD.  We  define  7  =  p0ff/pon  to 
describe  the  RDD  effect,  which  can  be  modeled  by  a  normal 
distribution  as  7'  =  p1  •  (1  +  a1  •  D).  Here  p1  =  Rh /Rl  and 
D  ~  1).  The  actual  value  of  o1  wifi  be  determined  by  the 

particular  device  structure,  material  and  fabrication  process.  In 
our  following  simulations,  we  assume  a1  =  2%. 

To  obtain  the  distributions  of  R'H  and  R'L,  we  conducted 
Monte-Carlo  simulations  with  10,000  3D  device  samples.  The 
variation  samples  are  generated  by  the  same  algorithm  and 

TABLE  I 

The  Ideal  Parameters  of  a  T1O2  Thin-film  Memristor. 


Geometry  Dimensions 

Electrical  Parameters 

Length  (/) 

50  nm 

Rh 

16,000  n 

Width  (w) 

50  nm 

Rl 

100  Q 

Thickness  ( h ) 

50  nm 

Uv 

10“14  m2S~1V~1 

Fig.  2.  The  Distribution  of  R'L  from  10,000  Monte-Carlo  simulations  and 
the  fitting  curve  of  Eq.  (11). 

constraints  in  [12].  Our  results  show  that  the  distributions  of 
both  R'h  and  R'L  are  close  to  normal  distributions  and  can  be 
approximated  by 

R'l  =  Rl  •  (rrl  +  ctRl  •  E),  (11) 

and 

R'h  =  Rh  '  (VRh  +  gRh  '  R)  ‘  (1  +  Gi  *  D)-  (12) 

Here,  two  independent  random  numbers  E  ~  A/*(0, 1)  and 
D  ~  ff(0, 1)  are  introduced.  E  represents  the  correlation 
between  R'H  and  R'L  due  to  the  same  geometry  variation 
sources.  D  represents  the  impact  of  RDD,  which  affects  the 
ratio  between  pQff  and  pon. 

Fig.  2  compares  the  approximated  normal  distribution 
shown  in  Eq.  (11)  and  the  actual  distribution  of  R'L  from 
Monte-Carlo  simulations.  The  root  mean  square  error  (RMSE) 
incurred  by  the  normal  distribution  approximation  is  only 
4.4%.  A  quantile-quantile  plot  is  also  shown  in  the  inset  of 
the  figure. 

B.  Distribution  of  a' 

As  shown  in  Eq.  (10),  the  calculation  of  the  average  doping 
front  a!  requires  time-consuming  filament-based  simulations. 
If  we  can  somehow  extract  the  actual  a'  directly  from  the 
designed  value  a  by  using  a  simplified  process  variation 
model,  the  simulation  cost  can  be  improved. 

Considering  the  fact  that  Rh  Rl  in  a  Ti02  memristor, 
a  simple  approximation  of  Eq.  (4)  can  be: 

q it)  =  1  -  VT^x,  (13) 

where, 

x  =  XjL  '(Jt  v(t)-dt  +  a 0  -  i  •  al).  (14) 

Eq.  (14)  shows  that  the  variation  of  a  can  be  directly 
linked  to  the  variation  of  X ,  which  has  three  independent 
contributors:  (1)  the  variation  of  thin-film  thickness  h ,  (2)  the 
impact  of  RDD  that  is  represented  by  7,  and  (3)  the  magnetic 
flux  of  the  input  signal  ip  =  fVdt.  Interestingly,  LER  does 
not  impact  a',  as  also  proved  by  the  simulations  in  [12]. 

Impact  of  flux  p  and  boundary  conditions. 
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»  3D,  -3sigma  ■  3D,  +3sigma  -Stat, -3sigma  -B-  Stat, +3sigma 


t(s) 

Fig.  4.  a'(t)  when  Vap  =  1. IV  based  on  100  Monte-Carlo  simulations. 

If  there  are  no  process  variations,  the  average  doping  front  a 
will  be  uniquely  determined  by  the  magnetic  flux  ip,  as  shown 
in  Eq.  (14).  However,  when  process  variations  are  considered, 
the  historical  profile  of  the  electrical  excitations  (instead  of 
only  the  absolute  value  of  ip)  will  introduce  the  additional 
variations  of  ip  by  interacting  with  the  device  variations  such 
as  thin  film  thickness  hi  and  RDD  7'. 

To  understand  how  the  historical  profiles  of  ip,  e.g.,  the  am¬ 
plitude  and  the  time  duration  etc.,  affect  the  variation  of  a ,  we 
conducted  the  Monte-Carlo  simulations  with  10,000  3D  device 
samples  and  traced  the  doping  front  the  position  a'(t)  in  every 
device  samples.  A  sinusoid  input  signal  V  =  Vap  •  sin  (27 rf  •  t) 
with  a  fixed  frequency  /  =  0.5 Hz  is  applied.  The  +3cr  and 
— 3cr  variations  of  a{t  =  Is),  which  has  been  normalized 
against  the  designed  value  of  a(t  =  Is)  when  varying  Vap 
from  0.1U  to  1.2V ,  are  shown  as  the  curves  labeled  with 
“3D,  +3sigma”  and  “3D,  -3sigma”  in  Fig.  3,  respectively.  Here 
ip(t  =  Is)  =  Vap/ (7 rf),  which  is  proportional  to  Vap. 

The  simulation  results  show  that  the  3cr  variance  of  a'  is 
approximately  proportional  to  ip,  except  when  a'  is  close  to  1 
(or  the  total  memristance  is  close  to  R'L ).  A  significant  over¬ 
shoot  of  the  +3cr  variance  at  Vap  =  1.1V  is  observed.  Based 
on  the  mathematic  expression  in  Eq.  (13-14),  it  indicates  that 
a  has  a  wider  distribution  when  X  is  approaching  1. 

A  closer  look  of  this  boundary  situation  is  shown  in  Fig. 
4.  It  records  the  movement  history  of  the  doping  fronts  of 
100  device  samples  when  Vap  =  1.1  TA  The  RED  curve  is 
the  designed  behavior  while  the  100  BLUE  curves  come  from 
Monte-Carlo  simulations.  The  average  doping  fronts  of  some 
devices  with  large  process  variations  hit  the  device  boundary 
during  the  movements  and  cause  the  large  variance  of  a'  at  the 
+3cr  corner.  However,  large  a'  variance  at  the  —  3cr  corner  not 
observed  in  Fig.  3.  It  is  because  the  average  doping  front  of 
the  most  devices  at  the  —  3cr  corner  will  not  reach  the  device 
boundary  under  the  simulated  electrical  excitations. 


TABLE  II 

Statistical  model  of  Ti02  thin-film  memristor. 


Variation  Parameters:  hrh,  &rh,  Hrl,  aRLi  07 
Coefficients:  w\,  w 2,  £1,  £2 
Independent  Random  Numbers: 

D  ~  A/"(0,  1),  E  ~  JV( 0,  1),  G  ~  A/"(0,  1). 

Memristance  boundary: 

R'h  —  Rh  •  {hrh  +  aRH  ■  E)  ■  (1  ary  ■  D) 
R'l  —  Rl  ■  (^rl  +  <?rl)  ■  E 

Input  flux: 

/  r  vap  ■  dt  (a'  <  1) 

(fieff  =  S  Ut0 

{  h2(RH  +  RL)/(2^VRL)  ( a '  >  1) 

Doping  front  position: 

Oi  =  p  ■  a 

^  (l  +  ¥3-ei+V?-e2  -  (W1  •  E-\-W2‘  G))  ■  (1  +  cr^,  •  D) 

Memristance: 

M' (a)  —  R'l  •  a'  +  R'h  ■  (1  —  Oi') 


Further  increasing  the  flux  ip,  i.e.,  increasing  Vap  from  1.1V 
to  1.2U  makes  the  variance  of  a'  drop,  as  also  shown  in  Fig. 

3.  Under  this  scenario,  the  amplitude  of  ip  is  so  large  that 

the  average  doping  front  of  the  most  simulated  devices  have 
reached  the  device  boundary.  These  devices  have  a  constant 
a'  =  1  and  therefore,  will  not  contribute  to  the  variances 
of  average  doping  fronts.  A  new  concept  named  “effective 
flux,”  which  includes  the  constraints  of  device  boundary,  can 
be  expressed  as 

m  ff  =  {  fto  VaP  '  dt  ( a  <  (15) 

^  ff  I  h\RH  +  RL)/(2fivRL)  (a'  >  1). 

Impact  of  process  variations. 

The  complexity  of  Eq.  (13)  and  (14)  make  it  infeasible  to 
derive  a  simple  analytical  expression  of  the  variance  of  a' 
even  if  we  assume  both  h'  and  7'  follow  normal  distributions. 
However,  we  are  still  able  to  construct  a  polynomial-based 
a'  model  to  approximate  the  variations  of  the  average  doping 
front  in  the  memristor  device. 

By  running  extensive  numerical  simulations  under  various 
conditions,  we  found  the  actual  a'  can  be  modeled  as  the 
product  of  the  designed  value  a  and  a  coefficient  77  that 
represents  the  influence  of  process  variations  as: 

a'  =  77  •  a.  (16) 

Here  77  can  be  expressed  by  a  heuristic  formula  (partially  from 
the  first  order  Taylor  expansion  of  Eq.  (13)  and  (14))  as 

^  (lJr<P-£\Jrip-£2-{w\-E-\-W2-G))-(l-\-ary-D)  '  (17) 

Similar  to  the  definitions  we  used  in  Section  III-A,  D  and 
E  are  two  independent  random  numbers  that  represent  the 
impact  of  the  RDD  and  the  geometry  variations,  respectively. 
To  avoid  overestimating  the  impact  of  geometry  variations  on 
a ',  a  new  random  number  G  ~  A/*(0, 1)  is  introduced  to  offset 
the  impact  of  LER.  e\  and  £2  are  two  scalars  extracted  from 
the  actual  simulations.  The  coefficients  w\  and  W2  represent 
the  weights  of  E  and  G ,  where  w\  +  w\  =  1. 
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C.  Distribution  of  M' 

By  modeling  R'L ,  R'H  and  a'  with  Eq.  (11)  (12)  (16),  the 
total  memristance  M'  of  a  Ti02  memristor  can  be  simply 
calculated  by  Eq.  (2),  where  the  design  values  of  Rr/Rl  and 
a  should  be  replaced  by  the  actual  values  R'H ,  R'L  and  a'. 
Table  II  summarizes  the  main  steps  and  equations  included  in 
our  proposed  statistical  model  of  TiC>2  memristors. 

IV.  Model  Generation  Flow 

In  this  section,  we  describe  the  methodology  to  extract  the 
parameters  included  in  our  proposed  statistical  model  of  Ti02 
memristors.  Some  critical  implementation  considerations  are 
also  discussed. 

A.  Parameter  Extraction 

Based  on  the  discussion  in  Section  III,  we  proposed  a  four- 
step  extraction  methodology  to  obtain  the  parameters  required 
by  our  proposed  statistical  model. 

Step  1:  Model  input  characterization. 

In  the  electrical  testing  of  a  memristor  device,  only  three  pa¬ 
rameters  can  be  measured  directly:  R'H ,  R'L  and  I(t)(orV(t). 
The  measurement  of  I(t ),  which  is  a  time-varying  variable, 
requires  two  doping  front  a  movements  following  0  — )►  1  and 
1  — y  0  directions  under  certain  electrical  excitations:  as  we 
discussed  in  Section  III-B,  the  electrical  excitations  must  be 
carefully  controlled  to  prevent  the  doping  front  from  hitting 
the  boundaries.  The  absolute  value  of  the  flux  ip  can  be  simply 
calculated  as  the  integral  of  the  applied  voltage  over  time. 

Step  2:  Distribution  generations  for  R'H  and  R'L. 

The  PDF  of  R'L  can  be  modeled  as  a  normal  distribution 
as  shown  in  Eq.  (11).  The  mean  prl  and  the  variance  grl 
can  be  easily  extracted  from  the  statistical  data  measured 
(or  simulated  if  Monte-Carlo  simulation  is  used)  in  Step  1. 
Though  R'h  does  not  strictly  follow  a  normal  distribution  in 
Eq.  (12),  a  mean  of  / lrh  can  still  be  extracted.  Since  being 
exposed  to  the  same  geometry  variations,  grh  ~  grl  in  Eq. 
(12).  As  we  shall  show  in  Section  V-A,  modeling  R'H  as  a 
normal  distribution  approximation  introduces  very  marginal 
discrepancies  from  the  statistical  data  and  incurs  very  little 
impact  on  the  accuracy  of  a  modeling.  Finally,  the  variance 
of  g1  needs  to  be  provided  by  material-level  characterizations. 

Step  3:  Distribution  generation  for  a'. 

For  every  memristor  sample,  the  total  memristance  and  their 
average  doping  front  position  at  a  given  time  stamp  tj  can  be 
calculated  as: 

M'itj)  =  (18) 

Monte-Carlo  simulations  can  derive  the  distribution  of 
M'{tj )  and  <x'(tj)  and  the  mean  pa  and  the  variance  aa  of 
a' (tj).  Then  two  scalars  -  e i  and  £2  can  be  extracted  by  Eq. 
(16)  and  (17)  with  estimation  method. 

We  proposed  a  heuristic  method  to  generate  w\  and  W2: 
We  start  with  an  assumption  that  w\  =  1  and  W2  =  0. 
The  derived  mean  pm  and  variance  gm  of  total  memristance 
M  are  then  compared  to  the  corresponding  parameters  of 


the  measured  data  M',  say,  /Pm  and  ctm.  If  the  impact  of 
geometry  variations  is  overestimated,  we  will  have  gm  >  crM 
and  jiM  7^  We  then  adjust  the  values  of  w\  and  W2  to 
reduce  the  mismatch.  The  optimal  approximations  of  w\  and 
W2  can  be  obtained  when  the  parameter  mismatch  is  below  an 
acceptable  threshold. 

Step  4:  Model  verification  and  improvement. 

The  above  parameters  only  take  a  few  minutes  to  be 
extracted  when  sufficient  data  provided  and  sophisticated 
estimation  method  applied.  After  all  necessary  parameters  are 
extracted,  the  statistical  models  can  be  constructed.  The  static 
and  dynamic  behaviors  of  TiC>2  memristors  under  various 
electrical  excitations  can  be  simulated  by  our  model  without 
conducting  the  time-consuming  Monte-Carlo  simulations.  The 
accuracy  of  our  proposed  model  can  still  be  verified  by  com¬ 
paring  to  the  Monte-Carlo  simulation  results.  The  accuracy 
of  our  proposed  model  could  be  improved  by  optimizing  the 
iteration  termination  conditions  of  the  parameter  extraction 
processes. 

B.  Flow  Summary 

The  pseudocode  of  the  simulation  flow  of  the  memristor 
dynamic  behavior  is  summarized  below: 


Algorithm  1:  Statistical  Model  for  Computer  Simulation 


E  =  JV(0,  1),  G  =  JV(0,  1),  D  =  JV(0,  1) 

generated  at  the  beginning  of  simulation; 

c pin  +  1)  =  ip(n)  +  V(n  +  l)-dt 
if  abs(ip(n  +  1))  >  ipeff 

tp(n  +  1)  ==  sign(ip(n  +  1))  •  <peff 

end  if 

V(n+1)  = 
vel{n  +  1)  =  iiv  ■ 


+  l)-£1+<p(n  +  l)-e2-(^ 
R-L  V(n+1) 

w 


’1  ■  E  +  w2-G))-  (1  +  CTry  ■  D) 


M  (n) 

a(n  +  1)  =  a(n )  +  vel(n  +  1)  •  dt 
ol  (n  +  1)  =  rj(n  +  1)  •  a(n  +  1) 

Check  boundary  condition  of  ol  (n  +  1) 

If  a'  (n  +  1)  reaches  boundary 
(p(n  +  1)  =  0 
end  if 

M' (n  +  1)  =  R'l  ■  ol  (n  -f  1)  +  R'h  ■  (1  —  ol  (n  +  1)) 

return  M 


V.  Model  verification  and  design  implication 
A.  Model  Verification 

Because  of  the  lack  of  published  data  on  memristor  device 
variations,  we  use  Monte-Carlo  simulation  results  as  our 
baseline  to  validate  the  proposed  statistical  model.  Three  basic 
memristor  testing  data  -  R'H,  R'L  and  /(£),  are  generated  from 
the  simulations  on  the  3D  device  structure  examples  generated 


TABLE  III 

Statistical  model  parameters 


Variation  Parameters 

Coefficients 

Urh  0.994 

£  i  -0.028 

U>rl  0.994 

e2  0.072 

ctrh  2.16% 

w\  0.98 

ctrl  2.16% 

W2  0.2 

cr7  2% 

TABLE  IV 

Variance  between  3D  model  and  statistical  model. 


a  —  t 

M  -  t 

V  -  I 

+3er 

4.97% 

2.03% 

2.20% 

—3a 

5.05% 

2.12% 

1.99% 
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- 3D,  +3sigma  - 3D,  -3sigma  - Stat, +3sigma  - Stat, -3sigma  — Theoratical 


0  0.5  1  1.5  2 

t(s) 

Fig.  5.  Comparison  of  a'(t)  at  ±3cr  corners. 


- 3D,  +3sigma  - 3D, -3sigma  - Stat,  +3sigma  - Stat,  -3sigma  - Theoratical 


- 3D, +3sigma  - 3D,  -3sigma  - Stat,  +3sigma  - Stat, -3sigma  Theoratical 


Fig.  7.  Comparison  of  I  —  V  at  ±3cr  corners. 

by  the  technique  in  [12].  However,  there  is  nothing  preventing 
us  from  using  the  measurement  data  the  input  of  our  model 
generation.  The  parameters  of  the  generated  statistical  model 
are  summarized  in  TABLE  III,  which  are  used  in  our  following 
simulations. 

In  this  section,  the  effectiveness  of  our  proposed  statistical 
model  is  validated  by  comparing  our  simulation  results  of 
the  different  electrical  properties  of  a  TiC>2  memristor  with 
that  from  the  Monte-Carlo  simulations  based  on  3D  device 
structure  samples.  Every  Monte-Carlo  simulation  run  include 
10,000  samples. 

Test  1:  Fixed  Input  Signal.  A  sinusoid  input  signal  V  = 
Vap  •  sin(27r/-£)  with  Vap  =  IV  and  /  =  0.5 Hz  is 
applied  to  the  simulated  memristor  device.  Fig.  5,  6,  and 
7  show  the  simulation  results  of  and  I  —  V 

characteristics,  respectively.  The  results  of  our  statistical  model 
(labeled  “stat”)  are  very  close  to  the  results  of  the  Monte- 
Carlo  simulations  with  3D  device  samples  (labeled  “3D”). 
For  example,  compared  to  3D  Monte-Carlo  simulations,  our 
memristance  results  shows  only  ^2%  discrepancy  over  the 
simulated  time  range  at  the  ±3 a  corners.  Table  IV  summarizes 
the  variances  of  each  electrical  property  at  each  corner. 


— • — 3D, -3sigma  ■  3D,  +3sigma  —  <*-  Stat,  -3sigma  —  □-  Stat,  +3sigma 
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4- 

o 
<u 

g  0% 

.2 

■>  -10% 

-20% 

0  2  4  6  8  10 

/(Hz) 

Fig.  8.  ±3cr  variances  of  M  vs.  /  at  t  =  Aj. 
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Fig.  9.  ±3<t  variances  of  M  vs.  Vap  at  t  =  Is. 


Test  2:  Frequency  Dependency.  The  same  sinusoid  input  signal 
is  applied  with  Vap  =  2V  while  /  changes  from  2 Hz 
to  10 Hz.  Again,  our  statistical  model  demonstrates  a  good 
approximation  over  the  whole  range,  as  shown  in  Fig.  8. 

Test  3:  Impact  of  Flux.  In  this  test,  we  keep  the  input  as  a 
sinusoid  signal  with  a  fixed  /  =  0.5 Hz.  We  vary  Vap  from 
0.1V  to  1.2V  and  measure  a  and  M  at  t  =  Is.  For  an  ideal 
device  without  any  process  variations,  a  reaches  to  1  when 
Vap  =  1.26  V,  which  is  the  upper  bound  of  our  simulation.  The 
corresponding  a  and  M  at  different  V^’s  are  shown  in  Fig.  3 
and  Fig.  9,  respectively.  Results  show  that  our  statistical  model 
fits  the  Monte-Carlo  simulation  results  very  well  except  when 
VaP  is  ~  1.1V.  The  reason  has  been  explained  in  Section  III-B. 
In  memristor-based  analog  circuit  design,  we  recommend  not 
to  operate  the  memristors  in  this  region  where  the  doping 
fronts  in  some  memristors  may  hit  the  device  boundary  under 
the  process  variations. 

Runtime  Comparison.  The  proposed  statistical  model  of  Ti02 
memristors  significantly  improves  runtime  cost:  10,000  device- 
based  Monte-Carlo  simulations  [12]  took  ~  2  hours  in  a 
MATLAB  environment,  while  the  same  number  of  simulations 
only  took  ~  2  seconds  by  using  our  proposed  statistical  model. 

B.  Applications 

The  statistical  TiC>2  memristor  model  can  be  used  to  analyze 
the  impact  of  process  variations  on  memristor-based  circuit 
designs.  In  this  section,  we  use  an  op-amp  based  variable  gain 
amplifier  (VGA)  [13] [14]  as  the  example  to  demonstrate  the 
effectiveness  of  our  model. 

1)  An  op -amp  based  VGA  with  one  memristor:  Fig.  10 
shows  a  simple  VGA  design  which  includes  an  op-amp,  a 
resistor  and  a  memristor  [13].  Vbias  is  the  signal  to  control 
the  memristor  programming.  The  VGA  gain  (G),  which  is  the 
ratio  between  the  output  signal  Vout  and  the  input  signal 
Vjn,  is  mainly  determined  by  the  memristance  Rm  and  the 
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based  VGA  with  one  memristor. 


resistance  Rp  as: 


n _  Vqut  _  Rm 

^  ~  VIN  ~  Rf  • 


(19) 


In  our  simulation,  we  set  the  resistor  Rp  =  1  Kil.  The 
memristance  (Rm)  of  an  ideal  memristor  device  sweeps  from 
Rl  =  1000  to  Rh  =  16KQ.  The  initial  state  of  the  memristor 
is  Rl  ,  at  which  the  domain  wall  is  at  the  position  of  a  =  1.  A 
low-frequency  sinusoidal  signal  Vbias  =  0.2V -sin^Tr*  0.5 -t) 
is  applied  to  control  the  memristance  of  the  memristor.  Vbias 
was  carefully  selected  in  order  to  drive  the  doping  front 
a  over  the  full  range  [0, 1]  and  demonstrate  the  impact 
of  boundary  conditions. A  high-frequency  sinusoidal  signal 
Vin  =  0.1V  •  sin(27r  •  IK  •  t)  is  applied  to  measure  G.  100 
Monte-Carlo  simulations  are  run  based  on  our  statistical  Ti02 
memristor  model. 

Fig.  11  shows  the  VGA  gains  from  the  Monte-Carlo  sim¬ 
ulations.  During  the  first  half  period  t  =  0s  ^  Is,  Vbias  is 
positive.  G  grows  as  Rm  increases  until  the  domain  wall  hits 
the  boundary  (a  =  0).  In  contract,  a  negative  Vbias  results 
in  the  decreases  of  Rm  and  G.  As  expected,  the  variation  of 
G  agrees  with  the  distribution  of  memristance.  For  example, 
the  ±3cr  of  G  at  t  =  Is  is  ~  ±7%,  which  is  consistent  with 
±3cr  variance  of  the  R'L  distribution  shown  in  Fig.  2. 

An  interesting  phenomenon  observed  in  Fig.  11  is:  during 
the  two  periods  that  the  memristance  increases  and  decreases, 
the  variation  of  G  shows  the  different  trends.  For  instance, 
the  VGA  has  the  same  /iq  at  t  =  0.6s  and  t  =  1.25s. 
However,  <JG(t=i .25s)  is  much  less  than  crG(t=o.6s)-  The  ex¬ 
planation  is  the  follows:  when  programming  memristor,  the 
memristance  changing  rate  of  the  memristors  with  high  total 
memristance  is  faster  than  the  one  of  the  memristors  with  low 
total  memristance.  Such  mechanism  spreads  the  distribution  of 
instantaneous  memristance  M  when  programming  M  to  high 
and  suppresses  M  distribution  when  programming  M  to  low. 

2)  An  op-amp  based  VGA  with  two  memristors :  We  can 
replace  the  resistor  Rp  in  Fig.  10  with  another  memristor  to 
achieve  more  programming  capability  in  VGA  designs.  If  we 


Fig.  13.  Series  of  100  Monte-Carlo  simulations  of  the  gains  of  the  op-amp 
based  VGA  with  two  memristors. 

connect  the  two  memristors  in  the  same  direction,  say,  the 
Ti02  portion  of  one  memristor  is  connected  to  the  Ti02_a; 
portion  of  the  other  one,  the  changes  of  their  memristances 
will  follow  the  same  trend.  Ideally,  the  gain  of  the  VGA  (G) 
will  keeps  constant  since  the  ratio  of  the  two  memristances 
doesn’t  change.  To  obtain  a  wide  programmable  gain  range 
of  the  VGA,  two  memristors  are  connected  back  to  back  as 
shown  in  Fig.  12. 

We  conducted  Monte-Carlo  simulations  with  the  statistical 
Ti02  memristor  model  by  assuming  Rmi  and  Rm2  have 
the  same  design  parameters.  In  the  simulation,  we  use  the 
same  input  signal  Vin  as  the  previous  design,  and  change 
the  amplitude  of  the  control  signal  Vbias  to  IV  in  order  to 
excise  two  memristors.  The  G  changes  over  time  when  Rmi 
and  Rm2  are  fully  uncorrelated  is  demonstrated  in  Fig.  13.  The 
input  signal  Vin ,  the  control  signal  Vbias »  the  current  through 
the  memristors  Im,  and  the  output  signal  Vout  are  shown  in 
Fig.  14.  Compared  to  the  VGA  with  one  memristor,  the  range 
of  the  programmable  gain  G  of  the  VGA  with  two  memristors 
is  significantly  increased,  following  with  the  increase  in  G 
variation:  +3 crG(t=2s)  =  22.4%  and  —3 crG(t=2s)  =  15.5%. 

Considering  that  the  two  memristors  are  made  by  the  same 
process  at  the  same  time,  there  must  be  some  correlations 
between  Rmi  and  Rm2-  Fig.  15  shows  ±3 ac  by  adjusting  the 
correlation  coefficient  of  the  two  memristors  from  0  to  1.  As 


of  the  op-amp  based  VGA  with  two  memristors. 
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Fig.  15.  The  impact  of  device  correlation  on  the  variance  of  VGA  gain. 
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Fig.  16.  Model  with  error  coefficient  B. 

expected,  the  simulation  results  show  that  a  larger  correlation 
results  in  a  smaller  G  variation.  When  the  two  memristors  are 
fully  correlated  (means  they  are  identical),  we  still  observe  ~ 
±2.8%  variation.  This  small  variations  come  from  the  intrinsic 
total  memristance  variations  of  the  two  memristors  and  cannot 
be  canceled  to  each  other  at  circuit  design  level. 

C.  Model  error  tolerance  analysis 

Our  previous  analysis  was  build  on  the  polynomial  repre¬ 
sentation  or  Gaussian  distribution  assumption  on  the  physical 
parameters  of  Ti02  memristors.  However,  ignoring  higher 
order  approximation  in  our  model  may  incur  inaccuracy  when 
the  variation  is  large.  As  an  example,  ionic  drift  in  a  Ti02 
memristor  has  an  exponential  relationship  with  the  applied 
electrical  field  [15].  Our  previous  linear  approximation  may 
not  accurately  describe  the  doping  front  movement  when  the 
applied  voltage  is  large. 

Based  on  Taylor  expansion,  a  higher  order  relationship  (e.g., 
the  second  order  term)  can  be  added  on  top  of  Eq.  (3).  The 
modified  velocity  of  doping  front  movement  v(t)  becomes 

«(*)  =  W=»v&-W&+B- .(j$|)2.  (20) 

Here  the  coefficient  B  stands  for  the  second  order  impact 
of  the  voltage/current  on  doping  front  movement,  which  is 
independent  to  other  process  variations.  The  normal  value  of 
B  is  between  [—0.1,  0.1]  [15].  We  take  the  exponential  model 
of  doping  front  movement  in  3D  Monte-Carlo  simulation 
and  replace  the  velocity  express  Eq.  (3)  with  Eq.  (20)  in 
our  statistical  model,  The  comparison  between  the  simulation 
results  of  both  approaches  for  the  variance  of  M  is  shown  in 
Fig.  16. 

As  we  can  see  from  the  simulation  results,  our  statistical 


model  companies  well  with  the  3D  device  model  when  B  is 
less  than  0.04.  When  B  is  larger  than  0.04,  the  results  start  to 
diverge  due  to  the  increased  weight  of  the  higher  order  term. 
Higher  order  approximation  will  help  alleviate  this  issue  at  the 
cost  of  model  complexity  and  run  time  increase. 

VI.  Conclusion 

In  this  work,  we  developed  a  statistical  model  for  Ti02 
thin-film  memristors  based  on  the  theoretical  understanding  of 
the  impact  of  process  variations.  The  corresponding  parameter 
extraction  strategy  and  implementation  considerations  are  also 
discussed.  Simulation  results  show  that  our  proposed  statis¬ 
tical  model  has  excellent  accuracy  and  significantly  reduced 
runtime  cost:  3^4  magnitudes  run  time  reduction  and  only 
~  2%  accuracy  degradation  are  achieved  compared  to  the 
existing  device-based  Monte-Carlo  simulations.  An  application 
of  memristor  based  op-amp  VGA  is  also  explored  with  our 
model. 
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